function SimulatedMoments  = v9_Moments(ParametersIn);

% Parameters to be estimated
gamma         =    ParametersIn(1);
alpha         =    ParametersIn(2);
c1            =    ParametersIn(3);
c2            =    ParametersIn(4);
b1            =    ParametersIn(5);
b2            =    ParametersIn(6);

z             =    3.6;

%Set parameters
p_a = 1;
p_b = 0;
p_c = 0.5;

% True reporting rates
q_a = 0.900;
q_b = 0.409;
q_c = 0.507;
q_d = 0.234;

%Perceived reporting rates (March 24)
q_hat_a = 0.88;
q_hat_b = 0.59;
q_hat_c = 0.77;
q_hat_d = 0.23;

%Truncated true reporting costs
k_a =  9.690 ;
k_b =  4.065 ;
k_c = 5.196 ;
k_d = 2.054 ;

%Truncated perceived reporting costs
k_hat_a = 9.461 ;
k_hat_b = 6.149 ;
k_hat_c =  8.205 ;
k_hat_d =  2.007  ;
 



m = 20 ;



% PS effort  if s=0

num2_a =  alpha*gamma*  ((q_hat_a*((m*p_a)-k_hat_a))+b1)  *  ((q_hat_a*m*(1-p_a))+b2);
num2_b =  alpha*gamma*  ((q_hat_b*((m*p_b)-k_hat_b))+b1)  *  ((q_hat_b*m*(1-p_b))+b2);
num2_c =  alpha*gamma*  ((q_hat_c*((m*p_c)-k_hat_c))+b1)  *  ((q_hat_c*m*(1-p_c))+b2);
num2_d =  alpha*gamma*  ((q_hat_d*(-k_hat_d))+b1)  *  b2;


den_a = (2*c1*c2) - ((gamma^2)* ((q_hat_a*((m*p_a)-k_hat_a))+b1)  *  ((q_hat_a*m*(1-p_a))+b2));
den_b = (2*c1*c2) - ((gamma^2)* ((q_hat_b*((m*p_b)-k_hat_b))+b1)  *  ((q_hat_b*m*(1-p_b))+b2));
den_c = (2*c1*c2) - ((gamma^2)* ((q_hat_c*((m*p_c)-k_hat_c))+b1)  *  ((q_hat_c*m*(1-p_c))+b2));
den_d = (2*c1*c2) - ((gamma^2)* ((q_hat_d*(-k_hat_d))+b1)  * b2);


e2n_a = num2_a / den_a;
e2n_b = num2_b / den_b;
e2n_c = num2_c / den_c;
e2n_d = num2_d / den_d;


% CHW effort if s=0


num1_a = alpha*((q_a*((m*p_a)-k_a))+b1)*c2;
num1_b = alpha*((q_b*((m*p_b)-k_b))+b1)*c2;
num1_c = alpha*((q_c*((m*p_c)-k_c))+b1)*c2;
num1_d = alpha*((q_d*(-k_d))+b1)*c2;


e1n_a = num1_a / den_a;
e1n_b = num1_b / den_b;
e1n_c = num1_c / den_c;
e1n_d = num1_d / den_d;


% PS effort if s>0

eta_a = (b1*z) +b2 + ( q_hat_a*((m* (1+ (p_a*(z-1)))-(k_hat_a*z) )) );
eta_b = (b1*z) +b2 + ( q_hat_b*((m* (1+ (p_b*(z-1)))-(k_hat_b*z) )) );
eta_c = (b1*z) +b2 + ( q_hat_c*((m* (1+ (p_c*(z-1)))-(k_hat_c*z) )) );
eta_d = (b1*z) +b2 + ( q_hat_d*( -(k_hat_d*z) ) );


dens_a = (8*z*c1*c2)-((gamma^2)*(eta_a^2));
dens_b = (8*z*c1*c2)-((gamma^2)*(eta_b^2));
dens_c = (8*z*c1*c2)-((gamma^2)*(eta_c^2));
dens_d = (8*z*c1*c2)-((gamma^2)*(eta_d^2));


nums2_a = gamma*alpha*(eta_a^2);
nums2_b = gamma*alpha*(eta_b^2);
nums2_c = gamma*alpha*(eta_c^2);
nums2_d = gamma*alpha*(eta_d^2);

e2s_a = nums2_a / dens_a;
e2s_b = nums2_b / dens_b;
e2s_c = nums2_c / dens_c;
e2s_d = nums2_d / dens_d;


% CHW effort if s>0

e1term1_a = (4*alpha*z*c2) / dens_a;
e1term1_b = (4*alpha*z*c2) / dens_b;
e1term1_c = (4*alpha*z*c2) / dens_c;
e1term1_d = (4*alpha*z*c2) / dens_d;


e1term2_a = (q_a*(eta_a + (2*q_hat_a*z*(k_hat_a-k_a))))/ (2*q_hat_a*z);
e1term2_b = (q_b*(eta_b + (2*q_hat_b*z*(k_hat_b-k_b))))/ (2*q_hat_b*z);
e1term2_c = (q_c*(eta_c + (2*q_hat_c*z*(k_hat_c-k_c))))/ (2*q_hat_c*z);
e1term2_d = (q_d*(eta_d + (2*q_hat_d*z*(k_hat_d-k_d))))/ (2*q_hat_d*z);

e1term3_a = ( b1 * (q_hat_a-q_a)) / q_hat_a ;
e1term3_b = ( b1 * (q_hat_b-q_b)) / q_hat_b ;
e1term3_c = ( b1 * (q_hat_c-q_c)) / q_hat_c ;
e1term3_d = ( b1 * (q_hat_d-q_d)) / q_hat_d ;

e1s_a = e1term1_a * (e1term2_a + e1term3_a);
e1s_b = e1term1_b * (e1term2_b + e1term3_b);
e1s_c = e1term1_c * (e1term2_c + e1term3_c);
e1s_d = e1term1_d * (e1term2_d + e1term3_d);







% Calculate condition for s=1 
condition_a = (b2+(q_hat_a*m)-(z*b1)+(q_hat_a*k_hat_a*z)) / (q_hat_a*m*(z+1));
condition_b = (b2+(q_hat_b*m)-(z*b1)+(q_hat_b*k_hat_b*z)) / (q_hat_b*m*(z+1));
condition_c = (b2+(q_hat_c*m)-(z*b1)+(q_hat_c*k_hat_c*z)) / (q_hat_c*m*(z+1));


% Calculate output
e1_d = e1n_d;
e2_d = e2n_d;
x_d = alpha* e1_d + (gamma * e1_d * e2_d);


if p_a <= condition_a
    e1_a = e1s_a;
    e2_a = e2s_a;
    x_a = alpha* e1_a + (gamma * e1_a * e2_a);
else
    e1_a = e1n_a;
    e2_a = e2n_a;
    x_a = alpha* e1_a + (gamma * e1_a * e2_a);

end


if p_b <= condition_b
    e1_b = e1s_b;
    e2_b = e2s_b;
    x_b = alpha* e1_b + (gamma * e1_b * e2_b);

else
    e1_b = e1n_b;
    e2_b = e2n_b;
    x_b = alpha* e1_b + (gamma * e1_b * e2_b);

end

if p_c <= condition_c
    e1_c = e1s_c;
    e2_c = e2s_c;
    x_c = alpha* e1_c + (gamma * e1_c * e2_c);

else
    e1_c = e1n_c;
    e2_c = e2n_c;
    x_c = alpha* e1_c + (gamma * e1_c * e2_c);

end




% Side payment
st_a = ((b2 + (q_hat_a*m) - (z*b1) - (q_hat_a*( (m*p_a*(1+z)) - (k_hat_a*z)) ) )) / (2*q_hat_a*z);
st_b = ((b2 + (q_hat_b*m) - (z*b1) - (q_hat_b*( (m*p_b*(1+z)) - (k_hat_b*z)) ) )) / (2*q_hat_b*z);
st_c = ((b2 + (q_hat_c*m) - (z*b1) - (q_hat_c*( (m*p_c*(1+z)) - (k_hat_c*z)) ) )) / (2*q_hat_c*z);


if p_a <= condition_a
    s_a = st_a;
else
   s_a = 0;
end

if p_b <= condition_b
    s_b = st_b;
else
   s_b = 0;
end

if p_c <= condition_c
    s_c = st_c;
else
   s_c = 0;
end

s_d =0;



t_a = s_a * x_a ;
t_b = s_b * x_b ;
t_c = s_c * x_c;

de2_a = e2_a - e2_d;
de2_b = e2_b - e2_d;
de2_c = e2_c - e2_d;

dx_a = x_a - x_d;
dx_b = x_b - x_d;
dx_c = x_c - x_d;

%SimulatedMoments    = [de2_a; de2_b; de2_c; e2_d; dx_a;  dx_b; dx_c; x_d ];
SimulatedMoments    = [e2_a; e2_b; e2_c; e2_d; x_a;  x_b; x_c; x_d ];

return
